** ---------------------------------- **
* Author: Abdulrazzak Tamim 																
* First version: 24/7/2020																
* Last modified: 26/4/2024							
** ---------------------------------- **

* User must define a global macro in order to run the program:					
* "path" contains the raw data																

clear all												
set more off												
program drop _all																	
set matsize 10000			

if "`c(username)'"=="abdulrazzaktamim" {															
gl path "/Users/abdulrazzaktamim/Library/CloudStorage/Dropbox/SoilDoc-Abdulrazzak"									
}

program main

	use "$path/replication/data/yields.dta", clear
	preserve
		drop if trt==5
		HH_behavior
		Adoption_By_Yield
		Normalize
		IPW
		ITT_IPW
		ITT
		ITT_Spec2
		HeterogeneityArea
		UQR_2016
		UQR_2016_IPW 
		HeterogeneityEC_IPW
		PoolingITT
		Attrition
	restore
	Normalize
	CrossVillageITT	

end


program HH_behavior
	
	g variable=""
	replace variable= "missing" if mi(q28_2019)
	replace variable= "no" if q28_2019==0 
	replace variable= "yes, maize" if q28_2019==1 & q29_2019==1
	replace variable= "yes, not maize" if q28_2019==1 & q29_2019==0
	ta variable
	ta variable q44_2019, m //Table 3//	
	
	
	count if variable == "missing" & mi(q44_2019)
	loc r1_1 = r(N)
	count if variable == "missing" & q44_2019=="didn't_cultivate"
	loc r1_2 = r(N)
	count if variable == "missinge" & q44_2019=="no"
	loc r1_3 = r(N)
	count if variable == "missing" & q44_2019=="yes"
	loc r1_4 = r(N)
	count if variable == "missing"
	loc r1_5 = r(N)	
	
	count if variable == "yes, maize" & mi(q44_2019)
	loc r2_1 = r(N)
	count if variable == "yes, maize" & q44_2019=="didn't_cultivate"
	loc r2_2 = r(N)
	count if variable == "yes, maize" & q44_2019=="no"
	loc r2_3 = r(N)
	count if variable == "yes, maize" & q44_2019=="yes"
	loc r2_4 = r(N)
	count if variable == "yes, maize"
	loc r2_5 = r(N)

	count if variable == "yes, not maize" & mi(q44_2019)
	loc r3_1 = r(N)
	count if variable == "yes, not maize" & q44_2019=="didn't_cultivate"
	loc r3_2 = r(N)
	count if variable == "yes, not maize" & q44_2019=="no"
	loc r3_3 = r(N)
	count if variable == "yes, not maize" & q44_2019=="yes"
	loc r3_4 = r(N)
	count if variable == "yes, not maize"
	loc r3_5 = r(N)

	count if variable == "no" & mi(q44_2019)
	loc r4_1 = r(N)
	count if variable == "no" & q44_2019=="didn't_cultivate"
	loc r4_2 = r(N)
	count if variable == "no" & q44_2019=="no"
	loc r4_3 = r(N)
	count if variable == "no" & q44_2019=="yes"
	loc r4_4 = r(N)
	count if variable == "no"
	loc r4_5 = r(N)
	
	count if mi(q44_2019)
	loc r5_1 = r(N)
	count if q44_2019=="didn't_cultivate"
	loc r5_2 = r(N)
	count if q44_2019=="no"
	loc r5_3 = r(N)
	count if q44_2019=="yes"
	loc r5_4 = r(N)
	count if !mi(variable)
	loc r5_5 = r(N)

	cap file close
	file open fh using "$path/replication/output/Table3.tex", replace write
	file write fh " Missing          & `r1_1' & `r1_2' & `r1_3' & `r1_4' & `r1_5' \\ " _n		
	file write fh " & Yes, maize     & `r2_1' & `r2_2' & `r2_3' & `r2_4' & `r2_5' \\ " _n		
	file write fh " & Yes, not maize & `r3_1' & `r3_2' & `r3_3' & `r3_4' & `r3_5' \\ " _n		
	file write fh " & No             & `r4_1' & `r4_2' & `r4_3' & `r4_4' & `r4_5' \\ " _n		
	file write fh " & Total          & `r5_1' & `r5_2' & `r5_3' & `r5_4' & `r5_5'    " _n		
	file close fh 	

	g same = q44_2019=="yes" if !mi(q44_2019)

	areg same i.trt2 if !mi(far_mz_kg_acre_q30_2019) & !mi(gps_2019), a(village) r cl(village)
	outreg2 using "$path/replication/output/TableA5.doc", replace bd(3) sd(3) nocons addtext(Village FE, YES) /// 
	label ctitle(cultivated same main maize plot (=1))
	
	erase "$path/replication/output/TableA5.txt"

end
 

program Adoption_By_Yield

	xtile quintile=far_mz_kg_acre_q30_2016, n(5)
	foreach x of varlist fert_dum_2* { 
		mat `x'X=J(5,3,.)
		g `x'_100 = `x'*100
		tabstat `x'_100 if fert_dum_2014!=. &fert_dum_2016!=. &fert_dum_2019!=. , by(quintile) stat(mean sd N) save
		return list
		mat `x'1 = r(Stat1)
		mat `x'2 = r(Stat2)
		mat `x'3 = r(Stat3)
		mat `x'4 = r(Stat4)
		mat `x'5 = r(Stat5)
		mat `x'T = r(StatTotal)
		mat list `x'T
		mat `x'X = `x'1' \ `x'2' \ `x'3' \ `x'4' \ `x'5' \ `x'T'
		matmap `x'X `x'N, map(round(@,0.01))
		mat list `x'N
		mat rownames `x'N = First Second Third Fourth Fifth All
		esttab matrix(`x'N) using "$path/replication/output/Table7.csv", ///
			   nonumber append cells(Mean(fmt(2)) SD(fmt(2)))
		drop `x'_100
	}

end


program Normalize

	***Normalized SR yields
	su far_mz_kg_acre_q30_2014 if trt2==1 & gps_2016!=.
	g norm_SR_q30_14 = (far_mz_kg_acre_q30_2014 -r(mean))/r(sd)
	su far_mz_kg_acre_q30_2016 if trt2==1 & gps_2016!=.
	g norm_SR_q30_16 = (far_mz_kg_acre_q30_2016 -r(mean))/r(sd)
	su far_mz_kg_acre_q30_2019 if trt2==1 & gps_2019!=.
	g norm_SR_q30_19 = (far_mz_kg_acre_q30_2019 -r(mean))/r(sd)

	***Normalized satellite
	su gps_2016 if trt2==1 & far_mz_kg_acre_q30_2016!=.
	g sat_norm_q30_16 = (gps_2016 -r(mean))/r(sd)
	su gps_2019 if trt2==1 & far_mz_kg_acre_q30_2019!=.
	g sat_norm_q30_19 = (gps_2019 -r(mean))/r(sd)
	replace sat_norm_q30_16=. if cultivate_2016==0
	replace sat_norm_q30_19=. if cultivate_2019==0

	***Normalized GPS yields
	su gps_mz_kg_acre_q30_2014 if trt2==1 & gps_2016!=.
	g norm_GPS_q30_14 = (gps_mz_kg_acre_q30_2014 -r(mean))/r(sd)
	su gps_mz_kg_acre_q30_2016 if trt2==1 & gps_2016!=.
	gen norm_GPS_q30_16 = (gps_mz_kg_acre_q30_2016 -r(mean))/r(sd)
	su gps_mz_kg_acre_q30_2019 if trt2==1 & gps_2019!=.
	g norm_GPS_q30_19 = (gps_mz_kg_acre_q30_2019 -r(mean))/r(sd)
	
		       ***********************
	********** keep relevant variables **********
	           ***********************
	keep trt2 trt vil respondent_id obj*_binary_201* mz_kg_q* far_mz_kg_acre_q* /// 
	    q44_2019 norm* gps_201* sat_norm* not_inter* gps_area* trt2_* /// 
		gps_mz_kg_acre_q* ecdscm_rec asstfullsc fert_dum_2* cultivate* /// 
		headage malehead headeduc headeducp hhsize asstfullsc_qrt
		
	su asstfullsc, det /*2014 assets*/
	g asset_p50=(asstfullsc>r(p50)) if asstfullsc!=.
	g good_EC = ecdscm_rec == 3 | ecdscm_rec==4 if ecdscm_rec!=.
	gl het good_EC

end


program IPW 

	g missing_2016 = mi(norm_SR_q30_16) | mi(gps_2016) 
	g missing_2019 = mi(norm_SR_q30_19) | mi(gps_2019)

	ta missing_2016, m 
	ta missing_2019, m 
	
	forv i=6(3)9 { 
		probit missing_201`i' i.trt##(c.headage i.malehead i.headeducp c.hhsize i.asstfullsc_qrt) i.village, r cl(village) asis
		outreg2 using "$path/replication/output/TableA6.doc", append bd(2) sd(2) nocons addtext(Village FE, YES) /// 
		label ctitle(sr_st_`i') /// 
		keep(2.trt##c.headage 2.trt##1.malehead 2.trt##1.headeducp 2.trt##c.hhsize 2.trt##2.asstfullsc_qrt 2.trt##3.asstfullsc_qrt 2.trt##4.asstfullsc_qrt 3.trt##c.headage 3.trt##1.malehead 3.trt##1.headeducp 3.trt##c.hhsize 3.trt##2.asstfullsc_qrt 3.trt##3.asstfullsc_qrt 3.trt##4.asstfullsc_qrt 4.trt##c.headage 4.trt##1.malehead 4.trt##1.headeducp 4.trt##c.hhsize 4.trt##2.asstfullsc_qrt 4.trt##3.asstfullsc_qrt 4.trt##4.asstfullsc_qrt)
		predict probhat_201`i'
		as inrange(probhat_201`i',0,1) if missing_201`i'!=.
		gen ipw_201`i' = 1/probhat_201`i'
	}
	
end


program ITT_IPW 

	mat A_6 = J(5,1,.)
	mat A_9 = J(5,1,.)
	mat B_6 = J(5,1,.)
	mat B_9 = J(5,1,.)
	mat C_6 = J(5,1,.)
	mat C_9 = J(5,1,.)
	mat D_6 = J(5,1,.)
	mat D_9 = J(5,1,.)
	 
	 areg norm_SR_q30_16 i.trt2 if gps_2016!=. [pw=ipw_2016], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
		su norm_SR_q30_16 if e(sample) & trt2==1
		mat A_6[1,1] = r(mean)
		mat A_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat A_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat A_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat A_6[5,1] = r(p)
		mat list A_6

	 areg sat_norm_q30_16 i.trt2 if far_mz_kg_acre_q30_2016!=. [pw=ipw_2016], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_16 if e(sample) & trt2==1
		mat B_6[1,1] = r(mean)
		mat B_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat B_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat B_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat B_6[5,1] = r(p)
		mat list B_6

	 areg norm_SR_q30_19 i.trt2 if gps_2019!=. [pw=ipw_2019], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_19 if e(sample) & trt2==1
		mat A_9[1,1] = r(mean)
		mat A_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat A_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat A_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat A_9[5,1] = r(p)
		mat list A_9

	 areg sat_norm_q30_19 i.trt2 if far_mz_kg_acre_q30_2019!=. [pw=ipw_2019], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_19 if e(sample) & trt2==1
		mat B_9[1,1] = r(mean)
		mat B_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat B_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat B_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat B_9[5,1] = r(p)
		mat list B_9
																			
	areg norm_SR_q30_16 i.trt2 i.not_inter_2016 i.obj_binary_2016 if gps_2016!=. [pw=ipw_2016], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_16 if e(sample) & trt2==1
		mat C_6[1,1] = r(mean)
		mat C_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat C_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat C_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat C_6[5,1] = r(p)
		mat list C_6

	areg sat_norm_q30_16 i.trt2 i.not_inter_2016 i.obj_binary_2016 if ///
	far_mz_kg_acre_q30_2016!=. [pw=ipw_2016], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_16 if e(sample) & trt2==1
		mat D_6[1,1] = r(mean)
		mat D_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat D_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat D_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat D_6[5,1] = r(p)
		mat list D_6

	areg norm_SR_q30_19 i.trt2 i.not_inter_2019 i.obj_binary_2019 if gps_2019!=. [pw=ipw_2019], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_19 if e(sample) & trt2==1
		mat C_9[1,1] = r(mean)
		mat C_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat C_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat C_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat C_9[5,1] = r(p)
		mat list C_9

	 areg sat_norm_q30_19 i.trt2 i.not_inter_2019 i.obj_binary_2019 if ///
	 far_mz_kg_acre_q30_2019!=. [pw=ipw_2019], a(village) r cl(village)
		outreg2 using "$path/replication/output/Table5.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_19 if e(sample) & trt2==1
		mat D_9[1,1] = r(mean)
		mat D_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat D_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat D_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat D_9[5,1] = r(p)
		mat list D_9
	erase "$path/replication/output/Table5.txt"
	
	*Export
	matmap A_6 A_6N, map(round(@,0.01))
	matmap A_9 A_9N, map(round(@,0.01))
	matmap B_6 B_6N, map(round(@,0.01))
	matmap B_9 B_9N, map(round(@,0.01))
	matmap C_6 C_6N, map(round(@,0.01))
	matmap C_9 C_9N, map(round(@,0.01))
	matmap D_6 D_6N, map(round(@,0.01))
	matmap D_9 D_9N, map(round(@,0.01))
	
	putexcel set "$path/replication/output/Table5_F_tests.xls", sheet("F_tests") replace
	putexcel B2=matrix(A_6N)
	putexcel C2=matrix(B_6N)
	putexcel D2=matrix(A_9N)
	putexcel E2=matrix(B_9N)
	putexcel F2=matrix(C_6N)
	putexcel G2=matrix(D_6N)
	putexcel H2=matrix(C_9N)
	putexcel I2=matrix(D_9N)

	putexcel B1 = "SR_16"
	putexcel C1 = "Sat_16"
	putexcel D1 = "SR_19"
	putexcel E1 = "Sat_19"
	
	putexcel F1 = "SR_16"
	putexcel G1 = "Sat_16"
	putexcel H1 = "SR_19"
	putexcel I1 = "Sat_19"
	
	putexcel A2 = "CONTROL MEAN"
	putexcel A3 = "CONTROL SD"
	putexcel A4 = "V_vs_R"
	putexcel A5 = "V_vs_RV"
	putexcel A6 = "R_vs_RV"
	
	putexcel A7 = "Controls"
	putexcel B7 = "NO"
	putexcel C7 = "NO"
	putexcel D7 = "NO"
	putexcel E7 = "NO"

	putexcel F7 = "YES"
	putexcel G7 = "YES"
	putexcel H7 = "YES"
	putexcel I7 = "YES"

end


program ITT

	 areg norm_SR_q30_16 i.trt2 if gps_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
		su norm_SR_q30_16 if e(sample) & trt2==1
		mat A_6[1,1] = r(mean)
		mat A_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat A_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat A_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat A_6[5,1] = r(p)
		mat list A_6

	 areg sat_norm_q30_16 i.trt2 if far_mz_kg_acre_q30_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_16 if e(sample) & trt2==1
		mat B_6[1,1] = r(mean)
		mat B_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat B_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat B_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat B_6[5,1] = r(p)
		mat list B_6

	 areg norm_SR_q30_19 i.trt2 if gps_2019!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_19 if e(sample) & trt2==1
		mat A_9[1,1] = r(mean)
		mat A_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat A_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat A_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat A_9[5,1] = r(p)
		mat list A_9

	 areg sat_norm_q30_19 i.trt2 if far_mz_kg_acre_q30_2019!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_19 if e(sample) & trt2==1
		mat B_9[1,1] = r(mean)
		mat B_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat B_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat B_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat B_9[5,1] = r(p)
		mat list B_9
																			
	areg norm_SR_q30_16 i.trt2 i.not_inter_2016 i.obj_binary_2016 if gps_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_16 if e(sample) & trt2==1
		mat C_6[1,1] = r(mean)
		mat C_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat C_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat C_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat C_6[5,1] = r(p)
		mat list C_6

	areg sat_norm_q30_16 i.trt2 i.not_inter_2016 i.obj_binary_2016 if ///
	far_mz_kg_acre_q30_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_16 if e(sample) & trt2==1
		mat D_6[1,1] = r(mean)
		mat D_6[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat D_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat D_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat D_6[5,1] = r(p)
		mat list D_6

	areg norm_SR_q30_19 i.trt2 i.not_inter_2019 i.obj_binary_2019 if gps_2019!=.,  a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_19 if e(sample) & trt2==1
		mat C_9[1,1] = r(mean)
		mat C_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat C_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat C_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat C_9[5,1] = r(p)
		mat list C_9

	 areg sat_norm_q30_19 i.trt2 i.not_inter_2019 i.obj_binary_2019 if ///
	 far_mz_kg_acre_q30_2019!=., a(village)	r cl(village)
		outreg2 using "$path/replication/output/TableA11.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_19 if e(sample) & trt2==1
		mat D_9[1,1] = r(mean)
		mat D_9[2,1] = r(sd)
		test 2.trt2=3.trt2
		mat D_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat D_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat D_9[5,1] = r(p)
		mat list D_9
	erase "$path/replication/output/TableA11.txt"
	
	*Export
	matmap A_6 A_6N, map(round(@,0.01))
	matmap A_9 A_9N, map(round(@,0.01))
	matmap B_6 B_6N, map(round(@,0.01))
	matmap B_9 B_9N, map(round(@,0.01))
	matmap C_6 C_6N, map(round(@,0.01))
	matmap C_9 C_9N, map(round(@,0.01))
	matmap D_6 D_6N, map(round(@,0.01))
	matmap D_9 D_9N, map(round(@,0.01))
	
	putexcel set "$path/replication/output/TableA11_F_tests.xls", sheet("F-tests") replace
	putexcel B2=matrix(A_6N)
	putexcel C2=matrix(B_6N)
	putexcel D2=matrix(A_9N)
	putexcel E2=matrix(B_9N)
	putexcel F2=matrix(C_6N)
	putexcel G2=matrix(D_6N)
	putexcel H2=matrix(C_9N)
	putexcel I2=matrix(D_9N)
	
	putexcel B1 = "SR_16"
	putexcel C1 = "Sat_16"
	putexcel D1 = "SR_19"
	putexcel E1 = "Sat_19"
	
	putexcel F1 = "SR_16"
	putexcel G1 = "Sat_16"
	putexcel H1 = "SR_19"
	putexcel I1 = "Sat_19"
	
	putexcel A2 = "CONTROL MEAN"
	putexcel A3 = "CONTROL SD"
	putexcel A4 = "V_vs_R"
	putexcel A5 = "V_vs_RV"
	putexcel A6 = "R_vs_RV"
	
	putexcel A7 = "Controls"
	putexcel B7 = "NO"
	putexcel C7 = "NO"
	putexcel D7 = "NO"
	putexcel E7 = "NO"

	putexcel F7 = "YES"
	putexcel G7 = "YES"
	putexcel H7 = "YES"
	putexcel I7 = "YES"

end


program ITT_Spec2

	g V=trt2==2|trt2==4 if trt2!=.
	g R=trt2==3|trt2==4 if trt2!=.

	 areg norm_SR_q30_16 i.V##i.R if gps_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
		keep(1.V 1.R 1.V#1.R)
		su norm_SR_q30_16 if e(sample) & trt2==1
		mat A_6[1,1] = r(mean)
		mat A_6[2,1] = r(sd)
		test 1.V=1.R
		mat A_6[3,1] = r(p)
		test 1.V=1.R#1.V
		mat A_6[4,1] = r(p)
		test 1.R=1.R#1.V
		mat A_6[5,1] = r(p)
		mat list A_6

	 areg sat_norm_q30_16 i.V##i.R if far_mz_kg_acre_q30_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)  /// 
		keep(1.V 1.R 1.V#1.R)
	 su sat_norm_q30_16 if e(sample) & trt2==1
		mat B_6[1,1] = r(mean)
		mat B_6[2,1] = r(sd)
		test 1.V=1.R
		mat B_6[3,1] = r(p)
		test 1.V=1.R#1.V
		mat B_6[4,1] = r(p)
		test 1.R=1.R#1.V
		mat B_6[5,1] = r(p)
		mat list B_6

	 areg norm_SR_q30_19 i.V##i.R if gps_2019!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
		keep(1.V 1.R 1.V#1.R)
	 su norm_SR_q30_19 if e(sample) & trt2==1
		mat A_9[1,1] = r(mean)
		mat A_9[2,1] = r(sd)
		test 1.V=1.R
		mat A_9[3,1] = r(p)
		test 1.V=1.R#1.V
		mat A_9[4,1] = r(p)
		test 1.R=1.R#1.V
		mat A_9[5,1] = r(p)
		mat list A_9

	 areg sat_norm_q30_19 i.V##i.R if far_mz_kg_acre_q30_2019!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
		keep(1.V 1.R 1.V#1.R)
	 su sat_norm_q30_19 if e(sample) & trt2==1
		mat B_9[1,1] = r(mean)
		mat B_9[2,1] = r(sd)
		test 1.V=1.R
		mat B_9[3,1] = r(p)
		test 1.V=1.R#1.V
		mat B_9[4,1] = r(p)
		test 1.R=1.R#1.V
		mat B_9[5,1] = r(p)
		mat list B_9
																			
	areg norm_SR_q30_16 i.V##i.R i.not_inter_2016 i.obj_binary_2016 if gps_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
		keep(1.V 1.R 1.V#1.R 1.not_inter_2016 1.obj_binary_2016)
	 su norm_SR_q30_16 if e(sample) & trt2==1
		mat C_6[1,1] = r(mean)
		mat C_6[2,1] = r(sd)
		test 1.V=1.R
		mat C_6[3,1] = r(p)
		test 1.V=1.R#1.V
		mat C_6[4,1] = r(p)
		test 1.R=1.R#1.V
		mat C_6[5,1] = r(p)
		mat list C_6

	areg sat_norm_q30_16 i.V##i.R i.not_inter_2016 i.obj_binary_2016 if ///
	far_mz_kg_acre_q30_2016!=., a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
		keep(1.V 1.R 1.V#1.R 1.not_inter_2016 1.obj_binary_2016)
	 su sat_norm_q30_16 if e(sample) & trt2==1
		mat D_6[1,1] = r(mean)
		mat D_6[2,1] = r(sd)
		test 1.V=1.R
		mat D_6[3,1] = r(p)
		test 1.V=1.R#1.V
		mat D_6[4,1] = r(p)
		test 1.R=1.R#1.V
		mat D_6[5,1] = r(p)
		mat list D_6

	areg norm_SR_q30_19 i.V##i.R i.not_inter_2019 i.obj_binary_2019 if gps_2019!=.,  a(village) r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
		keep(1.V 1.R 1.V#1.R 1.not_inter_2019 1.obj_binary_2019)
	 su norm_SR_q30_19 if e(sample) & trt2==1
		mat C_9[1,1] = r(mean)
		mat C_9[2,1] = r(sd)
		test 1.V=1.R
		mat C_9[3,1] = r(p)
		test 1.V=1.R#1.V
		mat C_9[4,1] = r(p)
		test 1.R=1.R#1.V
		mat C_9[5,1] = r(p)
		mat list C_9

	 areg sat_norm_q30_19 i.V##i.R i.not_inter_2019 i.obj_binary_2019 if ///
	 far_mz_kg_acre_q30_2019!=., a(village)	r cl(village)
		outreg2 using "$path/replication/output/TableA14.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
		keep(1.V 1.R 1.V#1.R 1.not_inter_2019 1.obj_binary_2019)
	 su sat_norm_q30_19 if e(sample) & trt2==1
		mat D_9[1,1] = r(mean)
		mat D_9[2,1] = r(sd)
		test 1.V=1.R
		mat D_9[3,1] = r(p)
		test 1.V=1.R#1.V
		mat D_9[4,1] = r(p)
		test 1.R=1.R#1.V
		mat D_9[5,1] = r(p)
		mat list D_9
	erase "$path/replication/output/TableA14.txt"
	
	*Export
	matmap A_6 A_6N, map(round(@,0.01))
	matmap A_9 A_9N, map(round(@,0.01))
	matmap B_6 B_6N, map(round(@,0.01))
	matmap B_9 B_9N, map(round(@,0.01))
	matmap C_6 C_6N, map(round(@,0.01))
	matmap C_9 C_9N, map(round(@,0.01))
	matmap D_6 D_6N, map(round(@,0.01))
	matmap D_9 D_9N, map(round(@,0.01))
	
	putexcel set "$path/replication/output/TableA14_F_tests.xls", sheet("F_tests") replace
	putexcel B2=matrix(A_6N)
	putexcel C2=matrix(B_6N)
	putexcel D2=matrix(A_9N)
	putexcel E2=matrix(B_9N)
	putexcel F2=matrix(C_6N)
	putexcel G2=matrix(D_6N)
	putexcel H2=matrix(C_9N)
	putexcel I2=matrix(D_9N)
	
	putexcel B1 = "SR_16"
	putexcel C1 = "Sat_16"
	putexcel D1 = "SR_19"
	putexcel E1 = "Sat_19"
	
	putexcel F1 = "SR_16"
	putexcel G1 = "Sat_16"
	putexcel H1 = "SR_19"
	putexcel I1 = "Sat_19"
	
	putexcel A2 = "CONTROL MEAN"
	putexcel A3 = "CONTROL SD"
	putexcel A4 = "V_vs_R"
	putexcel A5 = "V_vs_RV"
	putexcel A6 = "R_vs_RV"
	
	putexcel A7 = "Controls"
	putexcel B7 = "NO"
	putexcel C7 = "NO"
	putexcel D7 = "NO"
	putexcel E7 = "NO"

	putexcel F7 = "YES"
	putexcel G7 = "YES"
	putexcel H7 = "YES"
	putexcel I7 = "YES"

end

	
program HeterogeneityArea
	
	loc j=0
	forv i=0(.25)2 { 
		loc j=`j'+1
		areg sat_norm_q30_16 i.trt2 if norm_SR_q30_16!=. & not_inter_2016==1 ///
		     & gps_area_old_2016>`i' &gps_area_old_2016!=., a(village) r cl(village)
			estimates store Sat_`j'

		su gps_area_old_2016 if e(sample), det
		areg norm_SR_q30_16 i.trt2 if sat_norm_q30_16!=. & not_inter_2016==1 /// 
		     & gps_area_old_2016>`i' &gps_area_old_2016!=., a(village) r cl(village)	
			estimates store SR_`j'
		di `j'
	}	
	coefplot(Sat_1, color(red) symbol(O) label(Sat (all)) ciopts(lcolor(red))) ///
			(SR_1, color(blue) symbol(D) label(SR (all)) ciopts(lcolor(blue))) ///
			(Sat_7, color(red) symbol(O) label(Sat (>1.5)) ciopts(lcolor(red))) ///
			(SR_7, color(blue) symbol(D) label(SR (>1.5)) ciopts(lcolor(blue))) ///
			(Sat_8, color(red) symbol(O) label(Sat (>1.75)) ciopts(lcolor(red))) ///
			(SR_8, color(blue) symbol(D) label(SR (>1.75)) ciopts(lcolor(blue))) ///
			(Sat_9, color(red) symbol(O) label(Sat (>2)) ciopts(lcolor(red))) ///
			(SR_9, color(blue) symbol(D) label(SR (>2)) ciopts(lcolor(blue))), ///
			 vertical  yline(0) xline(5) ci(95) ylabel(-1(.25)1, nogrid) ///
			 coeflabels(2.trt2 = "Voucher" 3.trt2 = "Recommendations" 4.trt2 = "Voucher+Recommendations") ///
			 ytitle(Productivity in standard deviations (2016)) drop(_cons) /*legend(off)*/
	gr export "$path/replication/output/Figure4a.pdf", as(pdf) replace
	
	//with weights
	loc j=0
	forv i=0(.25)2 { 
		loc j=`j'+1
		areg sat_norm_q30_16 i.trt2 if norm_SR_q30_16!=. & not_inter_2016==1 ///
		     & gps_area_old_2016>`i' &gps_area_old_2016!=. [pw=ipw_2016], a(village) r cl(village)
			estimates store Sat_`j'

		su gps_area_old_2016 if e(sample), det
		areg norm_SR_q30_16 i.trt2 if sat_norm_q30_16!=. & not_inter_2016==1 /// 
		     & gps_area_old_2016>`i' &gps_area_old_2016!=. [pw=ipw_2016], a(village) r cl(village)	
			estimates store SR_`j'
		di `j'
	}	
	coefplot(Sat_1, color(red) symbol(O) label(Sat (all)) ciopts(lcolor(red))) ///
			(SR_1, color(blue) symbol(D) label(SR (all)) ciopts(lcolor(blue))) ///
			(Sat_7, color(red) symbol(O) label(Sat (>1.5)) ciopts(lcolor(red))) ///
			(SR_7, color(blue) symbol(D) label(SR (>1.5)) ciopts(lcolor(blue))) ///
			(Sat_8, color(red) symbol(O) label(Sat (>1.75)) ciopts(lcolor(red))) ///
			(SR_8, color(blue) symbol(D) label(SR (>1.75)) ciopts(lcolor(blue))) ///
			(Sat_9, color(red) symbol(O) label(Sat (>2)) ciopts(lcolor(red))) ///
			(SR_9, color(blue) symbol(D) label(SR (>2)) ciopts(lcolor(blue))), ///
			 vertical  yline(0) xline(5) ci(95) ylabel(-1(.25)1, nogrid) ///
			 coeflabels(2.trt2 = "Voucher" 3.trt2 = "Recommendations" 4.trt2 = "Voucher+Recommendations") ///
			 ytitle(Productivity in standard deviations (2016)) drop(_cons) /*legend(off)*/
	gr export "$path/replication/output/Figure4b.pdf", as(pdf) replace
	
end


program UQR_2016

	/*unconditional quantile regression: by Firpo et al. (2009)
	  recentered influence function (RIF) regressions */
	
	//Self-reported
	preserve
	forv i=10(5)90 {
		rifhdreg norm_SR_q30_16 trt2_2-trt2_4 if sat_norm_q30_16!=., a(village) rif(q(`i'))
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]

	}
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.0f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((a) Self-Reported) ytitle("Productivity in 2016 (s.d.)") ///
		legend(off) legend(rows(1)) name(SR_2016, replace)

	restore
	//Satellite
	clear matrix
	preserve
	forv i=10(5)90 {
		rifhdreg sat_norm_q30_16 trt2_2-trt2_4 if norm_SR_q30_16!=., a(village) rif(q(`i'))
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]

	}
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.1f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((b) Satellite-Based) ytitle("Productivity in 2016 (s.d.)") ///
		legend(off) legend(rows(1))	name(Sat_2016, replace)
	
	gr combine SR_2016 Sat_2016, cols(1) ysize(20) xsize(15)
	gr export "$path/replication/output/Figure3a.pdf", replace

	restore	
	
	
	//////////////////////////////
	/////////SR vs. GPS
	//////////////////////////////
	
	//Self-reported
	clear matrix
	preserve
	forv i=10(5)90 {
		rifhdreg far_mz_kg_acre_q30_2016 trt2_2-trt2_4 far_mz_kg_acre_q30_2014 if gps_mz_kg_acre_q30_2016!=. & gps_mz_kg_acre_q30_2014!=., a(village) rif(q(`i'))
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]

	}
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.0f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((a) Self-Reported) ytitle("Productivity in 2016 (kg/acre)") ///
		legend(off) legend(rows(1)) name(SR2_2016, replace)

	restore

	//GPS
	clear matrix
	preserve
	forv i=10(5)90 {
		rifhdreg gps_mz_kg_acre_q30_2016 trt2_2-trt2_4 gps_mz_kg_acre_q30_2014 if far_mz_kg_acre_q30_2016!=. & far_mz_kg_acre_q30_2014!=., a(village) rif(q(`i'))	
		//xtrifreg gps_mz_kg_acre_q30_2016 trt2_2-trt2_4 gps_mz_kg_acre_q30_2014 if far_mz_kg_acre_q30_2016!=. & far_mz_kg_acre_q30_2014!=., fe i(village) q(`i')		
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]
	}
	
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.0f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((b) GPS-Corrected) ytitle("Productivity in 2016 (kg/acre)") ///
		legend(off) legend(rows(1))	name(GPS_2016, replace)
	gr combine SR2_2016 GPS_2016, cols(1) ysize(20) xsize(15)
	gr export "$path/replication/output/Figure3b.pdf", replace
	restore	
	
end


program UQR_2016_IPW
	
	//Self-reported
	clear matrix
	preserve
	forv i=10(5)90 {
		rifhdreg norm_SR_q30_16 trt2_2-trt2_4 [pw=ipw_2016] if sat_norm_q30_16!=., a(village) rif(q(`i'))
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]

	}
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.0f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((a) Self-Reported) ytitle("Productivity in 2016 (s.d.)") ///
		legend(off) legend(rows(1)) name(SR_2016, replace)

	restore
	
	//Satellite
	clear matrix
	preserve
	forv i=10(5)90 {
		rifhdreg sat_norm_q30_16 trt2_2-trt2_4 [pw=ipw_2016] if norm_SR_q30_16!=., a(village) rif(q(`i'))
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]

	}
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.1f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((b) Satellite-Based) ytitle("Productivity in 2016 (s.d.)") ///
		legend(off) legend(rows(1))	name(Sat_2016, replace)
	
	gr combine SR_2016 Sat_2016, cols(1) ysize(20) xsize(15)
	gr export "$path/replication/output/Figure3c.pdf", replace

	restore	
	
	
	//////////////////////////////
	/////////SR vs. GPS
	//////////////////////////////
	
	//Self-reported
	clear matrix
	preserve
	forv i=10(5)90 {
		rifhdreg far_mz_kg_acre_q30_2016 trt2_2-trt2_4 far_mz_kg_acre_q30_2014 [pw=ipw_2016] if gps_mz_kg_acre_q30_2016!=. & gps_mz_kg_acre_q30_2014!=., a(village) rif(q(`i'))
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]

	}
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.0f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((a) Self-Reported) ytitle("Productivity in 2016 (kg/acre)") ///
		legend(off) legend(rows(1)) name(SR2_2016, replace)

	restore

	//GPS
	clear matrix
	preserve
	forv i=10(5)90 {
		rifhdreg gps_mz_kg_acre_q30_2016 trt2_2-trt2_4 gps_mz_kg_acre_q30_2014 [pw=ipw_2016] if far_mz_kg_acre_q30_2016!=. & far_mz_kg_acre_q30_2014!=., a(village) rif(q(`i'))		
		scalar beta = _b[trt2_4]
		scalar sigma = _se[trt2_4]
		
		scalar cil = beta - (1.96*sigma)
		scalar ciu = beta + (1.96*sigma)
		matrix UQR = [nullmat(UQR) \ `i', beta, cil, ciu]
	}
	
	clear
	svmat UQR
	ren UQR1 qtile
	ren UQR2 treatment
	ren UQR3 treatment_cilo
	ren UQR4 treatment_cihi
		
	g y0=0
	g Percentile=qtile
	sort Percentile
	replace y0=. if Percentile<10 | Percentile>90
	format treatment treatment_cihi treatment_cilo %6.0f 

	tw (line treatment Percentile, lwidth(medthick) lcolor(gray lpattern(line))) ///
	   (line y0 Percentile, lcolor(black)) ///
	   (line treatment_cihi Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)) ///
	   (line treatment_cilo Percentile, lwidth(thin) lpattern(longdash_dot_dot) lcolor(navy)), ///
	    xla(10 20 30 40 50 60 70 80 90, nogrid) yla(, nogrid) ///
		title((b) GPS-Corrected) ytitle("Productivity in 2016 (kg/acre)") ///
		legend(off) legend(rows(1))	name(GPS_2016, replace)
	gr combine SR2_2016 GPS_2016, cols(1) ysize(20) xsize(15)
	gr export "$path/replication/output/Figure3d.pdf", replace
	restore	
	
end


program HeterogeneityEC_IPW

	mat A_6_30 = J(5,1,.)
	mat B_6_30 = J(5,1,.)
	mat C_6_30 = J(5,1,.)
	mat A_9_30 = J(5,1,.)
	mat B_9_30 = J(5,1,.)
	mat C_9_30 = J(5,1,.)
	mat D_6_30 = J(3,1,.)
	mat E_6_30 = J(3,1,.)
	mat F_6_30 = J(3,1,.)
	mat D_9_30 = J(3,1,.)
	mat E_9_30 = J(3,1,.)
	mat F_9_30 = J(3,1,.)
	
	foreach x in $het {
		forval i=6(3)9 {
			forval j=30/30 {
			***SR
			areg norm_SR_q`j'_1`i' i.trt2##c.`x' if sat_norm_q`j'_1`i'!=. ///
										 & norm_GPS_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			su norm_SR_q`j'_1`i' if e(sample) & trt2==1
				mat A_`i'_`j'[1,1] = r(mean)
				mat A_`i'_`j'[2,1] = r(sd)
				test 2.trt2=3.trt2
				mat A_`i'_`j'[3,1] = r(p)
				test 2.trt2=4.trt2
				mat A_`i'_`j'[4,1] = r(p)
				test 3.trt2=4.trt2
				mat A_`i'_`j'[5,1] = r(p)
				test 2.trt2 + 2.t#c.`x' = 3.trt2 + 3.t#c.`x'
				mat D_`i'_`j'[1,1] = r(p)
				test 2.trt2 + 2.t#c.`x' = 4.trt2 + 4.t#c.`x'
				mat D_`i'_`j'[2,1] = r(p)
				test 3.trt2 + 3.t#c.`x' = 4.trt2 + 4.t#c.`x'
				mat D_`i'_`j'[3,1] = r(p)

			 outreg2 using "$path/replication/output/Table9_A.doc", ///
				append ctitle(SR1`i') bd(3) sd(3) nocons label addtext(Village FE, YES) /// 
				keep(2.trt2 3.trt2 4.trt2)
				
			*Marginal
			margins, expression( _b[2.t] + _b[2.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(SR1`i'V) bd(3) sd(3)
			qui areg norm_SR_q`j'_1`i' i.trt2##c.`x' if sat_norm_q`j'_1`i'!=. ///
										 & norm_GPS_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			margins, expression( _b[3.t] + _b[3.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(SR1`i'R) bd(3) sd(3)
			qui areg norm_SR_q`j'_1`i' i.trt2##c.`x' if sat_norm_q`j'_1`i'!=. ///
										 & norm_GPS_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			margins, expression( _b[4.t] + _b[4.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(SR1`i'RV) bd(3) sd(3)

			***Satellite
			qui areg sat_norm_q`j'_1`i' i.trt2##c.`x' if norm_SR_q`j'_1`i'!=. ///
										 & norm_GPS_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			su sat_norm_q`j'_1`i' if e(sample) & trt2==1
			mat B_`i'_`j'[1,1] = r(mean)
			mat B_`i'_`j'[2,1] = r(sd)
			test 2.trt2=3.trt2
			mat B_`i'_`j'[3,1] = r(p)
			test 2.trt2=4.trt2
			mat B_`i'_`j'[4,1] = r(p)
			test 3.trt2=4.trt2
			mat B_`i'_`j'[5,1] = r(p)

			test 2.trt2 + 2.t#c.`x' = 3.trt2 + 3.t#c.`x'
			mat E_`i'_`j'[1,1] = r(p)
			test 2.trt2 + 2.t#c.`x' = 4.trt2 + 4.t#c.`x'
			mat E_`i'_`j'[2,1] = r(p)
			test 3.trt2 + 3.t#c.`x' = 4.trt2 + 4.t#c.`x'
			mat E_`i'_`j'[3,1] = r(p)

			outreg2 using "$path/replication/output/Table9_A.doc", ///
				append ctitle(Sat1`i') nocons bd(3) sd(3) label addtext(Village FE, YES) /// 
				keep(2.trt2 3.trt2 4.trt2)
				
			*Marginal
			margins, expression( _b[2.t] + _b[2.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(Sat1`i'V) bd(3) sd(3)
			qui areg sat_norm_q`j'_1`i' i.trt2##c.`x' if norm_SR_q`j'_1`i'!=. ///
										 & norm_GPS_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			margins, expression( _b[3.t] + _b[3.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(Sat1`i'R) bd(3) sd(3)
			qui areg sat_norm_q`j'_1`i' i.trt2##c.`x' if norm_SR_q`j'_1`i'!=. ///
										 & norm_GPS_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			margins, expression( _b[4.t] + _b[4.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(Sat1`i'RV) bd(3) sd(3)

			***GPS
			qui areg norm_GPS_q`j'_1`i' i.trt2##c.`x' if sat_norm_q`j'_1`i'!=. ///
										 & norm_SR_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			 su norm_GPS_q`j'_1`i' if e(sample) & trt2==1
			 mat C_`i'_`j'[1,1] = r(mean)
			 mat C_`i'_`j'[2,1] = r(sd)
			 test 2.trt2=3.trt2
			 mat C_`i'_`j'[3,1] = r(p)
			 test 2.trt2=4.trt2
			 mat C_`i'_`j'[4,1] = r(p)
			 test 3.trt2=4.trt2
			 mat C_`i'_`j'[5,1] = r(p)

			 test 2.trt2 + 2.t#c.`x' = 3.trt2 + 3.t#c.`x'
			 mat F_`i'_`j'[1,1] = r(p)
			 test 2.trt2 + 2.t#c.`x' = 4.trt2 + 4.t#c.`x'
			 mat F_`i'_`j'[2,1] = r(p)
			 test 3.trt2 + 3.t#c.`x' = 4.trt2 + 4.t#c.`x'
			 mat F_`i'_`j'[3,1] = r(p)

			outreg2 using "$path/replication/output/Table9_A.doc", ///
				append ctitle(GPS1`i') nocons bd(3) sd(3) label addtext(Village FE, YES) /// 
				keep(2.trt2 3.trt2 4.trt2)
			*Marginal					
			margins, expression( _b[2.t] + _b[2.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(GPS1`i'V) bd(3) sd(3)
			qui areg norm_GPS_q`j'_1`i' i.trt2##c.`x' if sat_norm_q`j'_1`i'!=. ///
										 & norm_SR_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			margins, expression( _b[3.t] + _b[3.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(GPS1`i'R) bd(3) sd(3)
			qui areg norm_GPS_q`j'_1`i' i.trt2##c.`x' if sat_norm_q`j'_1`i'!=. ///
										 & norm_SR_q`j'_1`i'!=. [pw=ipw_201`i'], a(village) r cl(village)
			margins, expression( _b[4.t] + _b[4.t#c.`x']) post
			 outreg2 using "$path/replication/output/Table9_B.xls", ///
				append ctitle(GPS1`i'RV) bd(3) sd(3)
			}
		}
		
		*Export
		matmap A_6_30 A_6N_30, map(round(@,0.01))
		matmap B_6_30 B_6N_30, map(round(@,0.01))
		matmap C_6_30 C_6N_30, map(round(@,0.01))
		matmap A_9_30 A_9N_30, map(round(@,0.01))
		matmap B_9_30 B_9N_30, map(round(@,0.01))
		matmap C_9_30 C_9N_30, map(round(@,0.01))

		matmap D_6_30 D_6N_30, map(round(@,0.01))
		matmap E_6_30 E_6N_30, map(round(@,0.01))
		matmap F_6_30 F_6N_30, map(round(@,0.01))
		matmap D_9_30 D_9N_30, map(round(@,0.01))
		matmap E_9_30 E_9N_30, map(round(@,0.01))
		matmap F_9_30 F_9N_30, map(round(@,0.01))

		putexcel set "$path/replication/output/Table9_F_tests.xls", sheet("F-tests") replace
		putexcel B2=matrix(A_6N_30)
		putexcel C2=matrix(B_6N_30)
		putexcel D2=matrix(C_6N_30)
		putexcel E2=matrix(A_9N_30)
		putexcel F2=matrix(B_9N_30)
		putexcel G2=matrix(C_9N_30)
		putexcel B8=matrix(D_6N_30)
		putexcel C8=matrix(E_6N_30)
		putexcel D8=matrix(F_6N_30)
		putexcel E8=matrix(D_9N_30)
		putexcel F8=matrix(E_9N_30)
		putexcel G8=matrix(F_9N_30)
		
		putexcel B1 = "SR_16"
		putexcel C1 = "Sat_16"
		putexcel D1 = "GPS_16"
		putexcel E1 = "SR_19"
		putexcel F1 = "Sat_19"
		putexcel G1 = "GPS_19"
		putexcel A2 = "CONTROL MEAN"
		putexcel A3 = "CONTROL SD"
		putexcel A4 = "V_vs_R (Panel A)"
		putexcel A5 = "V_vs_RV (Panel A)"
		putexcel A6 = "R_vs_RV (Panel A)"
		putexcel A8 = "V_vs_R (Panel B)"
		putexcel A9 = "V_vs_RV (Panel B)"
		putexcel A10 = "R_vs_RV (Panel B)"
	
	}	
    erase "$path/replication/output/Table9_A.txt"
    erase "$path/replication/output/Table9_B.txt"
	
end


program PoolingITT

	preserve
		gen trt3=(trt2==2|trt2==4)
		
		mat A_6 = J(2,1,.)
		mat A_9 = J(2,1,.)
		mat B_6 = J(2,1,.)
		mat B_9 = J(2,1,.)
		mat C_6 = J(2,1,.)
		mat C_9 = J(2,1,.)
		mat D_6 = J(2,1,.)
		mat D_9 = J(2,1,.)
		mat E_6 = J(2,1,.)
		mat E_9 = J(2,1,.)
		mat F_6 = J(2,1,.)
		mat F_9 = J(2,1,.)
		mat G_6 = J(2,1,.)
		mat G_9 = J(2,1,.)
		mat H_6 = J(2,1,.)
		mat H_9 = J(2,1,.)

		 xtreg norm_SR_q30_16 i.trt3 if gps_2016!=., fe	r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
			su norm_SR_q30_16 if e(sample) & trt3==0
			mat A_6[1,1] = r(mean)
			mat A_6[2,1] = r(sd)
			mat list A_6

		 xtreg sat_norm_q30_16 i.trt3 if far_mz_kg_acre_q30_2016!=., fe r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
		 su sat_norm_q30_16 if e(sample) & trt3==0
			mat B_6[1,1] = r(mean)
			mat B_6[2,1] = r(sd)
			mat list B_6

		 xtreg norm_SR_q30_19 i.trt3 if gps_2019!=., fe	r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
		 su norm_SR_q30_19 if e(sample) & trt3==0
			mat A_9[1,1] = r(mean)
			mat A_9[2,1] = r(sd)
			mat list A_9

		 xtreg sat_norm_q30_19 i.trt3 if far_mz_kg_acre_q30_2019!=., fe	r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
		 su sat_norm_q30_19 if e(sample) & trt3==0
			mat B_9[1,1] = r(mean)
			mat B_9[2,1] = r(sd)
			mat list B_9

		 xtreg norm_SR_q30_16 i.trt3 i.not_inter_2016 i.obj_binary_2016 if gps_2016!=., ///
		 fe	r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
		 su norm_SR_q30_16 if e(sample) & trt3==0
			mat C_6[1,1] = r(mean)
			mat C_6[2,1] = r(sd)
			mat list C_6

		 xtreg sat_norm_q30_16 i.trt3 i.not_inter_2016 i.obj_binary_2016 if ///
		 far_mz_kg_acre_q30_2016!=., fe	r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
		 su sat_norm_q30_16 if e(sample) & trt3==0
			mat D_6[1,1] = r(mean)
			mat D_6[2,1] = r(sd)
			mat list D_6

		 xtreg norm_SR_q30_19 i.trt3 i.not_inter_2019 i.obj_binary_2019 if gps_2019!=., ///
		 fe	r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
		 su norm_SR_q30_19 if e(sample) & trt3==0
			mat C_9[1,1] = r(mean)
			mat C_9[2,1] = r(sd)
			mat list C_9

		 xtreg sat_norm_q30_19 i.trt3 i.not_inter_2019 i.obj_binary_2019 if ///
		 far_mz_kg_acre_q30_2019!=., fe	r cl(village)
			outreg2 using "$path/replication/output/TableA8.doc", ///
			append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
		 su sat_norm_q30_19 if e(sample) & trt3==0
			mat D_9[1,1] = r(mean)
			mat D_9[2,1] = r(sd)
			mat list D_9
		erase "$path/replication/output/TableA8.txt"
		
		*Export
		matmap A_6 A_6N, map(round(@,0.01))
		matmap A_9 A_9N, map(round(@,0.01))
		matmap B_6 B_6N, map(round(@,0.01))
		matmap B_9 B_9N, map(round(@,0.01))
		matmap C_6 C_6N, map(round(@,0.01))
		matmap C_9 C_9N, map(round(@,0.01))
		matmap D_6 D_6N, map(round(@,0.01))
		matmap D_9 D_9N, map(round(@,0.01))

		putexcel set "$path/replication/output/TableA8_F_tests.xls", sheet("F-tests") replace
		putexcel B2=matrix(A_6N)
		putexcel C2=matrix(B_6N)
		putexcel D2=matrix(A_9N)
		putexcel E2=matrix(B_9N)
		putexcel F2=matrix(C_6N)
		putexcel G2=matrix(D_6N)
		putexcel H2=matrix(C_9N)
		putexcel I2=matrix(D_9N)
		
		putexcel B1 = "SR_16"
		putexcel C1 = "Sat_16"
		putexcel D1 = "SR_19"
		putexcel E1 = "Sat_19"
		putexcel F1 = "SR_16"
		putexcel G1 = "Sat_16"
		putexcel H1 = "SR_19"
		putexcel I1 = "Sat_19"
		
		putexcel A2 = "CONTROL MEAN"
		putexcel A3 = "CONTROL SD"
		putexcel A4 = "Controls"
		putexcel B4 = "NO"
		putexcel C4 = "NO"
		putexcel D4 = "NO"
		putexcel E4 = "NO"
		putexcel F4 = "YES"
		putexcel G4 = "YES"
		putexcel H4 = "YES"
		putexcel I4 = "YES"
		
	restore
	
end


program Attrition

	mat T=J(1,6,.)
	
	forv i= 6(3)9 {
		forv j= 30/30 {
			g att_`j'_sr_201`i'  = (far_mz_kg_acre_q`j'_201`i'!=.)
			qui xtreg att_`j'_sr_201`i' i.trt2, fe r cl(village)
				outreg2 using "$path/replication/output/TableA3.doc", ///
				append ctitle(att_`j'_sr_201`i') bd(3) sd(3) nocons label addtext(Village FE, YES)
			su att_`j'_sr_201`i' if e(sample) & trt2==1
			loc control_sr_`i' = r(mean)
		
			g att_`j'_sat_201`i' = (sat_norm_q30_1`i'!=.) 
			qui xtreg att_`j'_sat_201`i' i.trt2, fe r cl(village)
				outreg2 using "$path/replication/output/TableA3.doc", ///
				append ctitle(att_`j'_sat_201`i') bd(3) sd(3) nocons label addtext(Village FE, YES)	
			su att_`j'_sat_201`i' if e(sample) & trt2==1
			loc control_sat_`i' = r(mean)
			
			g att_`j'_gps_201`i' = (gps_mz_kg_acre_q`j'_201`i'!=.)
			qui xtreg att_`j'_gps_201`i' i.trt2, fe r cl(village)
				outreg2 using "$path/replication/output/TableA3.doc", ///
				append ctitle(att_`j'_gps_201`i') bd(3) sd(3) nocons label addtext(Village FE, YES)	
			su att_`j'_gps_201`i' if e(sample) & trt2==1
			loc control_gps_`i' = r(mean)
			}
		}
		
	erase "$path/replication/output/TableA3.txt"

	mat T[1,1] = `control_sr_6'
	mat T[1,2] = `control_sat_6'
	mat T[1,3] = `control_gps_6'
	mat T[1,4] = `control_sr_9'
	mat T[1,5] = `control_sat_9'
	mat T[1,6] = `control_gps_9'
	matmap T T_new, map(round(@,0.01))
	mat list T_new 
	putexcel set "$path/replication/output/TableA3_mean.xls", sheet("Control") replace
	putexcel A2 = "control mean"
	putexcel B2 = matrix(T_new)

	mat A_30_6 = J(6,2,.)
	mat A_30_9 = J(6,2,.)
	mat S_30_6 = J(1,1,.)
	mat S_30_9 = J(1,1,.)
	mat G_30_6 = J(1,1,.)
	mat G_30_9 = J(1,1,.)
	
	loc k=0
	forval i= 6(3)9 {
	forval j= 30/30 {
		loc k = `k' + 1
		xtreg far_mz_kg_acre_q`j'_2014 att_`j'_sr_201`i'##trt2, fe r cl(village)
		outreg2 using "$path/replication/output/TableA4.doc", append ctitle ///
		(SR_`j') bd(2) sd(2) nocons label addtext(Village FE, YES) /// 
		keep(1.att_30_sr_201`i' att_30_sr_201`i'#2.trt2 2.trt2 /// 
		     1.att_30_sr_201`i' att_30_sr_201`i'#3.trt2 3.trt2 /// 
			 1.att_30_sr_201`i' att_30_sr_201`i'#4.trt2 4.trt2)
		su far_mz_kg_acre_q`j'_2014 if trt2==2 & att_`j'_sr_201`i'==0
		mat S_`j'_`i'[1,1]=r(mean)
		test (2.trt2) (1.att_`j'_sr_201`i'#2.trt2)
		mat A_`j'_`i'[`k',1] = r(p)
		mat list A_`j'_`i' 
		test (3.trt2) (1.att_`j'_sr_201`i'#3.trt2)
		loc k = `k' + 1
		mat A_`j'_`i'[`k',1] = r(p)
		mat list A_`j'_`i' 
		test (4.trt2) (1.att_`j'_sr_201`i'#4.trt2)
		loc k = `k' + 1
		mat A_`j'_`i'[`k',1] = r(p)
		mat list A_`j'_`i' 
		test (2.trt2) (1.att_`j'_sr_201`i'#2.trt2) (1.att_`j'_sr_201`i')
		loc k = `k' + 1
		mat A_`j'_`i'[`k',1] = r(p)
		mat list A_`j'_`i' 
		test (3.trt2) (1.att_`j'_sr_201`i'#3.trt2) (1.att_`j'_sr_201`i')
		loc k = `k' + 1
		mat A_`j'_`i'[`k',1] = r(p)
		mat list A_`j'_`i' 
		test (4.trt2) (1.att_`j'_sr_201`i'#4.trt2) (1.att_`j'_sr_201`i')
		loc k = `k' + 1
		mat A_`j'_`i'[`k',1] = r(p)
		mat list A_`j'_`i' 
	loc k=0
	xtreg gps_mz_kg_acre_q`j'_2014 att_`j'_gps_201`i'##trt2, fe r cl(village)
		su gps_mz_kg_acre_q`j'_2014 if trt2==2 & att_`j'_gps_201`i'==0
		mat G_`j'_`i'[1,1]=r(mean)
		outreg2 using "$path/replication/output/TableA4.doc", append ctitle ///
		(GPS_`j') bd(2) sd(2) nocons label addtext(Village FE, YES) /// 
		keep(1.att_30_gps_201`i' att_30_gps_201`i'#2.trt2 2.trt2 /// 
		     1.att_30_gps_201`i' att_30_gps_201`i'#3.trt2 3.trt2 /// 
			 1.att_30_gps_201`i' att_30_gps_201`i'#4.trt2 4.trt2)
		test (2.trt2) (1.att_`j'_gps_201`i'#2.trt2)
		loc k = `k' + 1
		mat A_`j'_`i'[`k',2] = r(p)
		mat list A_`j'_`i' 
		test (3.trt2) (1.att_`j'_gps_201`i'#3.trt2)
		loc k = `k' + 1
		mat A_`j'_`i'[`k',2] = r(p)
		mat list A_`j'_`i' 
		test (4.trt2) (1.att_`j'_gps_201`i'#4.trt2)
		loc k = `k' + 1
		mat A_`j'_`i'[`k',2] = r(p)
		mat list A_`j'_`i' 
		test (2.trt2) (1.att_`j'_gps_201`i'#2.trt2) (1.att_`j'_gps_201`i')
		loc k = `k' + 1
		mat A_`j'_`i'[`k',2] = r(p)
		mat list A_`j'_`i' 
		test (3.trt2) (1.att_`j'_gps_201`i'#3.trt2) (1.att_`j'_gps_201`i')
		loc k = `k' + 1
		mat A_`j'_`i'[`k',2] = r(p)
		mat list A_`j'_`i' 
		test (4.trt2) (1.att_`j'_gps_201`i'#4.trt2) (1.att_`j'_gps_201`i')
		loc k = `k' + 1
		mat A_`j'_`i'[`k',2] = r(p)
		mat list A_`j'_`i' 
		loc k=0
		}
	}

	mat T_30_6 = A_30_6
	mat T_30_9 = A_30_9
	matmap T_30_6 O_30_6, map(round(@,0.01))
	matmap T_30_9 O_30_9, map(round(@,0.01))
	putexcel set "$path/replication/output/TableA4_F_tests.xls", sheet("F-tests") replace
	putexcel B3=matrix(O_30_6)
	putexcel D3=matrix(O_30_9)
	mat A = S_30_6 \ G_30_6 \ S_30_9 \ G_30_9
	mat B = A'
	matmap B O1, map(round(@,0.01))
	putexcel B2=matrix(O1)
	
	putexcel B1 = "SR16"
	putexcel C1 = "GPS16"
	putexcel D1 = "SR19"
	putexcel E1 = "GPS19"
	putexcel A2 = "CONTROL MEAN"
	putexcel A3 = "IV-R: V"
	putexcel A4 = "IV-R: R"
	putexcel A5 = "IV-R: RV"
	putexcel A6 = "IV-P: V"
	putexcel A7 = "IV-P: R"
	putexcel A8 = "IV-P: RV"
	
	erase "$path/replication/output/TableA4.txt"
	
end


program CrossVillageITT

	mat A_6 = J(5,1,.)
	mat A_9 = J(5,1,.)
	mat B_6 = J(5,1,.)
	mat B_9 = J(5,1,.)
	mat C_6 = J(5,1,.)
	mat C_9 = J(5,1,.)
	mat D_6 = J(5,1,.)
	mat D_9 = J(5,1,.)

	 xtreg norm_SR_q30_16 i.trt2 if gps_2016!=., fe	r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	su norm_SR_q30_16 if e(sample) & trt2==1
	 	mat A_6[1,1] = r(mean)
		mat A_6[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat A_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat A_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat A_6[5,1] = r(p)
		mat list A_6

	 xtreg sat_norm_q30_16 i.trt2 if far_mz_kg_acre_q30_2016!=., fe r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_16 if e(sample) & trt2==1
		mat B_6[1,1] = r(mean)
		mat B_6[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat B_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat B_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat B_6[5,1] = r(p)
		mat list B_6

	 xtreg norm_SR_q30_19 i.trt2 if gps_2019!=., fe	r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_19 if e(sample) & trt2==1
		mat A_9[1,1] = r(mean)
		mat A_9[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat A_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat A_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat A_9[5,1] = r(p)
		mat list A_9

	 xtreg sat_norm_q30_19 i.trt2 if far_mz_kg_acre_q30_2019!=., fe	r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_19 if e(sample) & trt2==1
		mat B_9[1,1] = r(mean)
		mat B_9[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat B_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat B_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat B_9[5,1] = r(p)
		mat list B_9

	 xtreg norm_SR_q30_16 i.trt2 i.not_inter_2016 i.obj_binary_2016 if gps_2016!=., ///
	 fe	r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(SR16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_16 if e(sample) & trt2==1
	 	mat C_6[1,1] = r(mean)
		mat C_6[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat C_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat C_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat C_6[5,1] = r(p)
		mat list C_6

	 xtreg sat_norm_q30_16 i.trt2 i.not_inter_2016 i.obj_binary_2016 if ///
	 far_mz_kg_acre_q30_2016!=., fe	r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(Sat16) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_16 if e(sample) & trt2==1
	 	mat D_6[1,1] = r(mean)
		mat D_6[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat D_6[3,1] = r(p)
		test 2.trt2=4.trt2
		mat D_6[4,1] = r(p)
		test 3.trt2=4.trt2
		mat D_6[5,1] = r(p)
		mat list D_6

	 xtreg norm_SR_q30_19 i.trt2 i.not_inter_2019 i.obj_binary_2019 if gps_2019!=., ///
	 fe	r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(SR19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su norm_SR_q30_19 if e(sample) & trt2==1
	 	mat C_9[1,1] = r(mean)
		mat C_9[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat C_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat C_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat C_9[5,1] = r(p)
		mat list C_9

	 xtreg sat_norm_q30_19 i.trt2 i.not_inter_2019 i.obj_binary_2019 if ///
	 far_mz_kg_acre_q30_2019!=., fe	r cl(village)
	 	outreg2 using "$path/replication/output/TableA17.doc", ///
		append ctitle(Sat19) bd(3) sd(3) nocons label addtext(Village FE, YES)
	 su sat_norm_q30_19 if e(sample) & trt2==1
	 	mat D_9[1,1] = r(mean)
		mat D_9[2,1] = r(sd)
	 	test 2.trt2=3.trt2
		mat D_9[3,1] = r(p)
		test 2.trt2=4.trt2
		mat D_9[4,1] = r(p)
		test 3.trt2=4.trt2
		mat D_9[5,1] = r(p)
		mat list D_9
	erase "$path/replication/output/TableA17.txt"

	*Export
	matmap A_6 A_6N, map(round(@,0.01))
	matmap A_9 A_9N, map(round(@,0.01))
	matmap B_6 B_6N, map(round(@,0.01))
	matmap B_9 B_9N, map(round(@,0.01))
	matmap C_6 C_6N, map(round(@,0.01))
	matmap C_9 C_9N, map(round(@,0.01))
	matmap D_6 D_6N, map(round(@,0.01))
	matmap D_9 D_9N, map(round(@,0.01))

	putexcel set "$path/replication/output/TableA17_F_tests.xls", sheet("F-tests") replace
	putexcel B2=matrix(A_6N)
	putexcel C2=matrix(B_6N)
	putexcel D2=matrix(A_9N)
	putexcel E2=matrix(B_9N)
	putexcel F2=matrix(C_6N)
	putexcel G2=matrix(D_6N)
	putexcel H2=matrix(C_9N)
	putexcel I2=matrix(D_9N)
	
	putexcel B1 = "SR_16"
	putexcel C1 = "Sat_16"
	putexcel D1 = "SR_19"
	putexcel E1 = "Sat_19"
	
	putexcel F1 = "SR_16"
	putexcel G1 = "Sat_16"
	putexcel H1 = "SR_19"
	putexcel I1 = "Sat_19"
	
	putexcel A2 = "CONTROL MEAN"
	putexcel A3 = "CONTROL SD"
	putexcel A4 = "V_vs_R"
	putexcel A5 = "V_vs_RV"
	putexcel A6 = "R_vs_RV"
	
	putexcel A7 = "Controls"
	putexcel B7 = "NO"
	putexcel C7 = "NO"
	putexcel D7 = "NO"
	putexcel E7 = "NO"

	putexcel F7 = "YES"
	putexcel G7 = "YES"
	putexcel H7 = "YES"
	putexcel I7 = "YES"
	
	di 0.444 - 1.96*0.171
	di 0.444 + 1.96*0.171
	di 0.053 - 1.96*0.088
	di 0.053 + 1.96*0.088
	di 0.223 - 1.96*0.111
	di 0.223 + 1.96*0.111
	di 0.039 - 1.96*0.101
	di 0.039 + 1.96*0.101

end


*Execute													
main

*EOF
